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ABSTRACT 
(N 

^ The inner edge of the classical habitable zone is often defined by the critical flux needed to trigger the runaway greenhouse instability. 

This ID notion of a critical flux, however, may not be so relevant for inhomogeneously irradiated planets, or when the water content 
is limited (land planets). 

Here, based on results from our 3D global climate model, we present general features of the climate and large scale circulation on 
QQ close-in terrestrial planets. We find that the circulation pattern can shift from super-rotation to stellar/anti stellar circulation when 

the equatorial Rossby deformation radius significantly exceeds the planetary radius, changing the redistribution properties of the 
atmosphere. Using analytical and numerical arguments, we also demonstrate the presence of systematic biases between mean surface 
I— I temperatures or temperature profiles predicted from either ID or 3D simulations. 

Including a complete modeling of the water cycle, we further demonstrate that for land planets closer than the inner edge of the 
classical habitable zone, two stable climate regimes can exist. One is the classical runaway state where all the water is vaporized, and 
• the other is a collapsed state where water is captured in permanent cold traps. We identify this "moist" bistability as the result of a 

competition between the greenhouse effect of water vapor and its condensation on the night side or near the poles, highlighting the 
dynamical nature of the runaway greenhouse effect. We also present synthetic spectra showing the observable signature of these two 
JL states. 

^ Taking the exampfe of two prototype pfanets in this regime, namely Gl 581 c and HD 85512 b, we argue that depending on the rate of 

water delivery and atmospheric escape during the life of these planets, they could accumulate a significant amount of water ice at their 
^ surface. If such a thick ice cap is present, various physical mechanisms observed on Earth (e.g. gravity driven ice flows, geothermal 

1 I flux) should come into play to produce long-lived liquid water at the edge and/or bottom of the ice cap. Consequently, the habitability 

of planets at smaller orbital distance than the inner edge of the classical habitable zone cannot be ruled out. Transiting planets in this 
T-H regime represent promising targets for upcoming exoplanet characterization observatories like EChO and JWST. 

> 

l> 
O 

1. Beyond the runaway greenhouse limit equilibrium balance and is heated until the surface can radiate 

Pf-v at optical wavelength around a temperature of HOOK (Kasting 

y Because of the bias of planet characterization observatories to- ^934; Kasting 1988). All the water at the surface is va- 

^ ward strongly irradiated objects on short orbits, planets identified -^^^ ^^^^^ ^ greenhouse state. Roughly 

^ to be near the inner edge of the habitable zone represent valuable gpg^king, this radiation limit stems from the fact that the surface 

^ targets. This however raises numerous questions concerning our ^^^^ -^^^^^^^ -^^ temperature to radiate more, but the amount of 

^ ^ understanding of the processes delimiting this inner edge. ^^^^^ ^^p^^. troposphere, and thus the opacity of the latter, 

k> For a planet like the Earth, which from the point of view of increases exponentially with this suri^ace temperature (because 

climate can be considered as an aqua planet', a key constraint ^j^^ Clausius-Clapeyron relation). Then the thermal emission 

is provided by the positive feedback of water vapor on climate, ^f the opaque atmosphere is emitted at the level at which the op- 

When the insolation of the planet is increased, more water vapor ^j^al depth is close to unity, a level whose temperature reaches an 

IS released in the atmosphere by the hotter surface, increasing asymptotic value, hence the hmited flux^ (Nakajima et al. 1992; 

the greenhouse effect. If on Eaith today this effect is balanced pierrehumbert 2010) 

by other processes, in particular the increase of infrared emission ,^,.1 1 , • • j- j- ^■ 

, , J. • , 11 1 11 While large uncertainties remain regarding the radiative 

by the surface with temperature. It has been shown that above a . r . j fi, i- 

■ \c, .1 1 I - effect of water vapor continuum and the nonlinear, three- 
critical absorbed stellar flux, the planet cannot reach radiative ,. . irr^^j • jij.u -. ii c 

^ dimensional eftect of dynamics and clouds, the critical value of 
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' i.e. a planet whose surface is largely covered by connected oceans - A similar but slightly subtler limit stems from considering the radia- 

that have a planet-wide impact on the climate. Not to be confused with tive/thermodynamic equilibrium of the stratosphere (Kombayashi 1967; 

an ocean planet for which, in addition, water must represent a significant Ingersoll 1969; Nakajima et al. 1992), but it is less constraining in the 

fraction of the bulk mass (Leger et al. 2004). case considered. 
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the average absorbed stellar radiation, (1 - A)F*/4 (where 
is the stellar flux at the substellar point and A is the planetary 
bond albedo), should conservatively range between 240 W/m- 
(the mean absorbed stellar radiation on Earth) and 350 W/m^ 
(Abe et al. 2011). Terrestrial planets emitting such (low) IR 
fluxes will be difficult to characterize in the near future. 

Fortunately, the runaway greenhouse limit discussed above 
is not as fundamental as it may appear Indeed, the theory out- 
lined above assumes that a sufficient^ reservoir of liquid water is 
available everywhere on the surface (Selsis et al. 2007). While 
it would seem from ID modeling that liquid water would disap- 
pear at a lower surface temperature for a planet with a limited 
water inventory, it goes the opposite way in 3D. Global climate 
models show that such "land planets" can have a wider habitable 
zone (Abe et al. 201 1). Near the inner edge, this property stems 
from the fact that atmospheric transport of water from warm to 
cold regions of the planets is not counteracted by large-scale sur- 
face runoff' (Abe et al. 2005). Heavily irradiated regions are then 
drier and can emit more IR ffux than the limit found for a satu- 
rated atmosphere. 

In principle, it is thus possible to find planets absorbing and 
re-emitting more ffux than the runaway greenhouse threshold 
and for which liquid/solid water can be thermodynamically sta- 
ble at the surface, provided that they accreted limited water sup- 
plies (or lost most of them). As will be shown hereafter, this 
could be the case for G1581 c (Udry et al. 2007; Mayor et al. 
2009; Forveille et al. 2011) and HD 85512b (Pepe et al. 2011), 
the only two confirmed planets having masses compatible with 
a terrestrial planet (Mp < 7 - SMg; although defining such a 
mass limit is questionable), lying beyond the inner edge of the 
classical habitable zone but receiving a moderate stellar ffux, 
F^^./4 < 1000 W/m^. Note that these planets lie even beyond the 
inner edge of the limits set by Abe et al. (201 1) for land planets. 

However, for such close-in planets whose rotation states have 
been strongly affected by tidal dissipation and which are either 
in a state of synchronous, pseudo-synchronous'* or resonant rota- 
tion with a near zero obliquity (Goldreich & Peale 1966; Heller 
et al. 2011; Makarov et al. 2012), efficient and permanent cold 
traps present on the dark hemisphere or at the poles could irre- 
versibly capture all the available water in a permanent ice cap. 
To tackle this issue, we have performed several sets of global cli- 
mate model simulations for prototype land planets with the char- 
acteristic masses and (strong) incoming stellar ffuxes of G1581 c 
and HD 85512 b. 3D climate models are necessary to assess the 
habitability of land planets in such scenarios because horizontal 
inhomogeneities in impinging ffux and water distribution are the 
key features governing the climate. This may explain why this 
scenario has received little attention so far. 

One of the difficulties lies in the fact that, as will be exten- 
sively discussed here, there is no precise ffux limit above which 
a runaway greenhouse state is inevitable. Instead, for a wide 
range of parameters, two stable equilibrium climate regimes ex- 
ist because of the very strong positive feedback of water vapor. 
The actual regime is determined not only by the incoming stel- 
lar ffux and spectrum, or by the planet's atmosphere, but also by 
the amount of water available and its initial distribution. Con- 
sequently, instead of trying to cover the whole parameter space, 
here we have performed simulations for two already detected 



' Meaning that, if fully vaporized, the water reservoir must produce a 
surface pressure which is higher than the pressure of the critical point 
of water (220bars). 

Note that the stability of the pseudo-synchronous rotation state is still 
debated for terrestrial planets (Makarov & Efroimsky 2012) 



planets. We will show that qualitatively robust features appear 
in our simulations, and argue that the processes governing the 
climate discussed here are quite general. 

In Sect. 2 we introduce the various components of our nu- 
merical climate model and we present the basic climatology of 
our model in a completely dry case (Sect. 3). We show that the 
atmospheric circulation can settle into two different regimes - 
namely super-rotation or stellar/anti-stellar circulation - and that 
the transition of one regime to the other occurs when the equa- 
torial Rossby deformation radius significantly exceeds the plan- 
etary radius. We also show that strong discrepancies can occur 
between ID and globally averaged 3D results. Then we show 
that when water is present, the climate regime reached depends 
on the initial conditions and can exhibit a bistability (Sect. 4). 
We argue that this bistability is a dynamical manifestation of the 
interplay between the runaway greenhouse instability on the hot 
hemisphere and the cold trapping arising either on the dark side 
or near the poles. In Sect. 5, we discuss the long-term stability 
of the different regimes found here and show that, if atmospheric 
escape is efficient, it could provide a stabilizing feedback for the 
climate. In Sect. 6, we suggest that if the planet has been able to 
accumulate a large amount of ice in its cold traps such that ice 
draining by glaciers is significant, liquid water could be present 
at the edge of the ice cap where glaciers melt. We also dis- 
cuss the possible presence of subsurface water under a thick ice 
cap. Finally, in Sect. 7 we present synthetic spectra and discuss 
the observable signatures of the two climate regimes discussed 
above. 



2. Method 

Our simulations have been performed with an upgraded version 
of the LMD generic Global Climate Model (GCM) specifically 
developed for the study of extrasolar planets (Wordsworth et al. 
2010, 2011; Selsis et al. 2011) and paleoclimates (Wordsworth 
et al. 2013; Forget et al. 2013). The model uses the 3D dynam- 
ical core of the LMDZ 3 GCM (Hourdin et al. 2006), based on 
a finite-difference formulation of the primitive equations of geo- 
physical ffuid dynamics. A spatial resolution of 64 x 48 x 20 in 
longitude, latitude and altitude is used for the simulations. 

2.1. radiative transfer 

The method used to produce our radiative transfer model is sim- 
ilar to Wordsworth et al. (2011). For a given mixture of at- 
mospheric gases, here N2 with 376ppmv of CO2 and a vari- 
able amount of water vapor, we computed high resolution spec- 
tra over a range of temperatures and pressures using the HI- 
TRAN 2008 database (Rothman et al. 2009). Because the CO2 
concentration depends on the efficiency of complex weathering 
processes, it is mostly unconstrained for hot synchronous land 
planets. The choice of an Earth-like concentration is thus ar- 
bitrary but our prototype atmosphere should be representative 
of a background atmosphere with a low greenhouse effect. For 
this study we used temperature, pressure grids with values T = 
{110, 170, ...,710) K,/? = {10-^10-^..., lO^lmbar. The H2O 
volume mixing ratio could vary in the range {10^^, 10"^, 1). 
H2O lines were truncated at 25 cm while the water vapor con- 
tinuum was included using the CKD model (Clough et al. 1989). 
We also account for opacity due to N2-N2 collision induced ab- 
sorption (Borysow & Frommhold 1986; Richard et al. 201 1). 

The correlated-k method was then used to produce a smaller 
database of coefficients suitable for fast calculation in a GCM. 
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Planet name 




G1581C 


HD 85512 b 


Stellar luminosity 




0.0135 


0.126 


Stellar mass 




0.31 


0.69 


Orbital semi-major axis 


a [AU] 


0.073 


0.26 


Orbital rotation rate 


□orb [S-'] 


5.59x10-^ 


1.25x10-^ 


Orbital eccentricity 


e 


0-0.05 


0-0.11 


Obliquity 


ep 





- 


Tidal resonance 


n 


1, 3:2 


- 


atmospheric pressure 


p, [bar] 


0.1-10 


- 


Mass"'" 


Mp [Me] 


6.25 


4.15 


Radius* 




1.85 


1.60 


Surface gravity 


g [m s-^] 


18.4 


15.8 


Surface albedo 




0.3 




Ice albedo 




0.3-0.55 




Surface roughness 


zo [m] 


1 X 10-2 




Specific heat capacity 


Cp [JK-'kg 


-1] 1003.16 




Condensation nuclei 


A^c [kg-'] 


1 X 10^ 





' This is the mass used in simulations. It is inferred from the Mp sin / 
measured by radial velocity by assuming / = 60°. This inclination 
seems favored by recent observations of a disk around G1581. 
* Inferred from the mass assuming Earth density. 

Table 1. Standard parameters used in the climate simulations. 



Thanks to the linearity of the Schwarzschild equation of radia- 
tive transfer, the contribution of the thermal emission and down- 
welling stellar radiation can be treated separately, even in the 
same spectral channel. Hence, we do not assume any spec- 
tral separation between the stellar and planetary emission wave- 
lengths. For thermal emission, the model uses 19 spectral bands, 
and the two-stream equations are solved using the hemispheric 
mean approximation (Toon et al. 1989). Infrared scattering by 
clouds is also taken into account. Absorption and scattering of 
the downwelling stellar radiation is treated with the 5-Eddington 
approximation within 18 bands. Sixteen points are used for the 
^-space integration, where g is the cumulated distribution func- 
tion of the absorption data for each band. 

As in Wordsworth et al. (2010), the emission spectrum of our 
M dwarf, Gl 581, was modeled using the Virtual Planet Labora- 
tory AD Leo data (Segura et al. 2005). For the K star HD 85512, 
we used the synthetic spectrum of a K dwarf with an effective 
temperature of 4700 K, g - lO^-^ ms"^ ^nd [M/H] - -0.5 dex 
from the database of Allard et al. (2000). 

2.2. water cycle 

In the atmosphere, we follow the evolution of water in its vapor 
and condensed phase. These tracers are advected by the dynam- 
ical core, mixed by turbulence and dry and moist convection. 

Much care has been devoted to develop a robust and nu- 
merically stable water cycle scheme that is accurate both in the 
trace gas (water vapor mass mixing ratio q'v <k 1) and dom- 
inant gas (^v ~ 1) limit. In particular, the atmospheric mass 
and surface pressure variation (and the related vertical trans- 
port of tracers through pressure levels) due to any evapora- 
tion/sublimation/condensation of water vapor is taken into ac- 
count with a scheme similar to the one developed by Forget et al. 
(2006) for CO2 atmospheres. 

Cloud formation is treated using the prognostic equations of 
Le Treut & Li (1991). For each column and level, this scheme 



provides the cloud fraction, /c, and the mass mixing ratio amount 
of condensed water, q\, as a function of and the saturation 
vapor mass mixing ratio, ^s- In addition, when part of a col- 
umn reaches both 100% saturation and a super-moist-adiabatic 
lapse rate, moist convective adjustment is performed following 
Manabe & Wetherald (1967) and the cloud fraction is set to 
unity. In the cloudy part of the sky, the condensed water is 
assumed to form a number Nc of cloud particles per unit mass 
of moist air This quantity, which represents the number den- 
sity per unit mass of activated cloud condensation nuclei, is as- 
sumed to be uniform throughout the atmosphere (Forget et al. 
2013). Then, the effective radius of the cloud particles is given 

by ^eff - 1^1 ^jip^ , where pc is the density of condensed 

water (10^ kg/m^ for liquid and 920 kg/m^ for ice). Precipitation 
is computed with the scheme of Boucher et al. (1995). As this 
scheme explicitly takes into account the dependence on gravity, 
cloud particle radii and background air density, it should remain 
valid for a wide range of situations. Finally, the total cloud frac- 
tion of the column is assumed to be the cloud fraction of the 
level with the thickest cloud, and radiative transfer is computed 
both in the cloudy and clear sky regions. Fluxes are then linearly 
weighted between the two regions. 

On the surface, ice can also have a radiative effect by lin- 
early increasing the albedo of the ground to Ai^e (see Table 1 ; in 
the present version of the code, both surface and ice albedos are 
independent of the wavelength) until the ice surface density ex- 
ceeds a certain threshold (here 30 kg m-^). In the general case, 
the values taken for both of these parameters can have a sig- 
nificant impact on the mean radiative balance of a planet. In 
addition, providing an accurate ice mean albedo can prove diffi- 
cult, as it depends on the ice cover, light incidence, state of the 
ice/snow, but also on the stellar spectrum (ice absorbing more 
efficiently infrared light; Joshi & Haberle 2012). However, in 
our case of a strongly irradiated, spin synchronized land planet 
where water is in limited supply, most of the day side is above 
0°C, and ice can form only in regions receiving little or no vis- 
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ible flux. As a result, and as confirmed by our tests, the value 
chosen for Aice has a negligible impact on the results presented 
hereafter 



3. Climate on synchronous planets: what 1D cannot 
tell you 

So far, most studies concerned with the determination of the lim- 
its of the habitable zone have used ID calculations. While these 
single column models can provide reasonable answers for plan- 
ets with (very) dense atmospheres and/or a rapid rotation that 
limits large-scale temperature contrasts, we will show here that 
only a three-dimensional model can treat the case of synchro- 
nized or slowly rotating exoplanets properly. To that purpose we 
will quantify both analytically and numerically the discrepancies 
that can arise between ID and 3D models. 

3.1. Equilibrium temperatures vs piiysical temperatures 

When considering the energetic balance of an object, one of the 
first quantities coming into mind is the equilibrium temperature, 
i.e. the temperature that a blackbody would need to have to radi- 
ate the absorbed flux. If one assumes local thermal equilibrium 
of each atmospheric column, one can define the local equilib- 
rium temperature 



f(l -A)F^n^ 



1/4 



(1) 



by Eq. (2) . Rigorously, this stems from the fact that T^eq,^^ 
is a concave function of the absorbed flux which, with the 
help of the Rogers-Holder inequality (Rogers 1888; Holder 
1889), yields 



1/4 



eq- 



(3) 



In other words, regions receiving twice as much flux do not 
need to be twice as hot to reach equilibrium and cold regions 
have a more important weight in the mean temperature. 

To give an idea of how strong this second effect can be, let 
us consider the limiting case of an airless, synchronous planet in 
thermal equilibrium with a constant surface albedo, A. In this 
case, the local surface temperature is equal to the local equilib- 
rium temperature and the mean surface temperature is given by 



471 

4-n 



dyU* 



(l-A)n//* 



cr 



1/4 



1/4 ^2;r 



r d0 r ,T 

Jo Jo 



(4) 



Performing the integrals, and using Eq. (2), we find 



2 V2 

r« = —-T„ 



0.57 Te, 



(5) 



where is the stellar flux at the substellar point at the time 
considered, A is the local bond albedo, ju* is the cosine of the 
stellar zenith angle (assumed equal to zero when the star is below 
the horizon), and cr is the Stefan-Boltzmann constant. 

Relaxing the local equilibrium condition but still assuming 
global thermal equilibrium of the planet, one can define the ef- 
fective equilibrium temperature 



- eq 



Hi' 



{\-A)F^^i^AS 



1/4 



(l-A)F^ 
4o- 



1/4 



(2) 



where S is the planet's surface area and A is the effective, plan- 
etary bond albedo. 

It is important to note that, because Tgq is derived from global 
energetic considerations, it is not necessarily close to the mean 
surface temperature, T^, and is often not. This is for two main 
reasons: 

• The most widely recognized of these reasons is that, when an 
absorbing atmosphere is present, the emitted flux does not di- 
rectly come from the ground. Because of this greenhouse ef- 
fect, the surface temperature of an isolated atmospheric col- 
umn can be higher than its equilibrium temperature. This 
explains why Venus and Earth mean surface temperatures 
are 737 K and 288 K while their effective equilibrium tem- 
peratures are around 23 1 K and 255 K (assuming an albedo 
of 0.75 and 0.29). 

• The second reason is the non-linearity of the blackbody bolo- 
metric emission with temperature. Because this emissions 
scales as T'^, it can be demonstrated that the spatial average 
of the local equilibrium temperature (Eq. (1)) is necessarily 
smaller than the effective equilibrium temperature defined 



While any energy redistribution by the atmosphere will tend to 
lower this difference*", one needs to remember that even the mean 
surface temperature can be much lower than the effective equi- 
librium temperature. 

These simple considerations show that equilibrium temper- 
ature can be misleading and should be used with great care. 
Indeed, if ID models can take into account the greenhouse ef- 
fect, they cannot predict energy redistribution on the night side. 
Consequently, not only can they not predict temperatures on the 
nightside, but they also predict fairly inaccurate mean tempera- 
tures. Thus, while equilibrium temperature remains a handy tool 
to compare the temperature regime from one planet to another, 
it is not accurate enough to draw conclusions on the climate of a 
given object. 



' Similar conclusions can be reached concerning the temporal average 
of the equilibrium temperature, which is smaller than the equilibrium 
temperature computed for the temporally averaged flux. Whether this 
behavior has an important impact on the physical temperature, how- 
ever, depends on the ratio of the thermal equilibration timescale of the 
atmosphere and ground to the period of the irradiation variation. For 
the Earth, thanks to the oceans, the thermal inertia is large enough so 
that the mean flux seems to be the relevant quantity for both diurnal and 
seasonal variations. 

* Note that if energy redistribution afl"ects the mean surface tempera- 
ture, it has no impact on the mean equilibrium temperature. The 1/4 fac- 
tor found in Eq. (2) and the inequality in Eq. (3) are exact and do not rely 
on any assumption about this redistribution or the thickness of the atmo- 
sphere. In fact, Eq. (3) can be generalized to temperature means at any 

pressure levels in an atmosphere by } / T{p)AS < / T\p)AS) , 
although this does not have any implications for the fluxes at the given 
level. In particular the inequality fj < T^q has no reason to hold when 
a sufficiently opaque atmosphere is present because of the greenhouse 
effect; Venus being a perfect example. 
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Fig. 1. Surface temperature maps (in °C) for the dry Gl 581 c (top row) and HD 85512 b (bottom row) runs for a total pressure (pj of 0.1, 1 and 
10 bars (from left to right). Longitude 0° corresponds to substellar point. 



3.2. Climate on a dry synchronous planet 

We now present a first set of 3D simulations of G1581 c and 
HD 85512 b that are modeled as completely dry planets. The 
goal of these idealized simulations is to gain insight into the var- 
ious circulation patterns arising on synchronously rotating ter- 
restrial planets and how the circulation regime depends on the 
atmospheric mass (mean surface pressure) and planetary rota- 
tion rate. 



3.2.1. Surface temperatures 

These dry runs have been performed using the parameters from 
Table 1, without a water cycle. For various surface pressures, 
the model has been initialized from a rest state with uniform 
surface temperature and integrated until statistical equilibrium is 
reached. 

The temperature maps for three sample surface pressures for 
the two planets are shown in Fig. 1. The most apparent feature 
is the expected dichotomy between the sunlit (between -90° and 
90° longitude) and dark hemispheres. A surprising outcome of 
these simulations is the great magnitude of the temperature dif- 
ference between the substellar point and the night side compared 
to the one in synchronous Earth simulations of Joshi (2003) and 
Merlis & Schneider (2010) or to the equator-pole temperature 
difference on Earth. To understand these differences, we have 
performed an additional suite of simulations of fast rotating or 
tidally synchronized Earth with our full radiative transfer or with 
a gray gas approximation (as in Merlis & Schneider 2010) and 
with or without water cycle. The model recovers both Earth 
present climatology and results similar to Merlis & Schneider 
(2010). 

Thanks to this suite of simulations, the causes of this impor- 
tant temperature contrast are found to be (more or less by order 
of importance): 



1. Latent cooling/heating. When water is abundantly present, 
the temperature at the substellar point is controlled mostly 
thermodynamically and is much lower than in the dry case. 
Latent energy transported by the atmosphere to the night side 
also contributes significantly to the warming of the latter. 

2. Non gray opacities. In our case dry case, greenhouse is 
solely due to the I5fj.m band of CO2 (especially on the night 
side). Thus, using a gray opacity of the background atmo- 
sphere with an optical depth near unity, as done in Merlis & 
Schneider (2010), considerably underestimates the cooling 
of the surface and thus overestimates the night side tempera- 
ture. 

3. A different planet. When its bulk density is constant, both 
the size and the surface gravity increase with the planet's 
size. Apart from the dynamical effect due to the variation 
of the Rossby number (extensively discussed below), note 
that for a given surface pressure (/?s), both the mass of the 
atmosphere and its radiative timescale (xi-ad) decrease with 
gravity (and with incoming flux for the latter) and that the 
advective timescale (radv) increases with the radius. For a 
given surface pressure the energy transport, thought to scale 
with the atmospheric mass and Tiad/Tadv thus tends to be less 
efficient on a more massive and more irradiated object. 

Indeed, as expected, the temperature difference decreases 
with increasing total pressure, because the strength of the ad- 
vection of sensible heat to the night side increases with atmo- 
spheric mass. Both because of this more efficient redistribution 
and of the greenhouse effect of the atmosphere, the globally aver- 
aged mean surface temperature increases with surface pressure, 
as shown in Fig. 2. 

For comparison, we also computed the surface temperature 
yielded by the ID version of our model. In this case, the atmo- 
spheric column is illuminated with a constant flux equal to the 
globally averaged impinging flux (1/4 of the substellar value) 
and with a stellar zenith angle of 60°. Apart from this prescrip- 
tion for the flux and the fact that large-scale dynamics is turned 
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100 



0.1 0.2 0.5 1.0 2.0 5.0 10.0 

Ps (bar) 



Fig. 2. Globally averaged (solid, filled) and minimal (dashed, empty) 
temperatures as a function of surface pressure for G1581 c (black cir- 
cles) and HD 85512 b (gray squares). For comparison, the surface tem- 
perature of the ID model for G1581 c is also shown (dotted curve). The 
error bars on the 200 mbar case show the temperatures reached by vary- 
ing the surface albedo between 0.1 and 0.5. 



off, our ID code is exactly the same as our 3D one. As can 
be seen in Fig. 2, even with a lObar atmosphere, the ID model 
is still 30 K hotter than the globally averaged 3D temperature. 
This difference increases up to 70 K for the 100 mb case because 
redistribution is less efficient (see Sect. 3.1). Note that a surface 
albedo of 0.3 is assumed here, similar to that in desert regions on 
Earth. The sensitivity of the surface temperature when albedo is 
varied between 0.1 and 0.5 is shown in Fig. 2. 

A troubling, and quite counterintuitive outcome of these 
simulations is that, as shown in Fig. 2, although G1581c re- 
ceives more flux and has a higher mean surface temperature than 
HD 85512 b, its minimal temperatures are always lower As can 
be seen in Fig. 1, the models for G1581 c present very cold re- 
gions located at mid latitudes and west from the western termina- 
tor (also called morning limb), whereas the HD 85512 b models 
have a rather uniform (and higher in average) temperature over 
the whole nightside. 

3.2.2. Circulation regime 

Further investigations show that this difference is a signature of 
the distinct circulation regimes arising on the two planets. In- 
deed, as shown in Fig. 3, our simulations of the atmosphere of 
G1581 c develop eastward jets (super-rotation) as found in many 
GCM of tidally locked exoplanets (Showman et al. 2008; Thras- 
tarson & Cho 2010; Heng & Vogt 201 1; Heng et al. 2011; Show- 
man et al. 2011; Showman & Polvani 2011; Showman et al. 
2013). Conversely, the atmosphere of HD 85512 b seems to pre- 
dominantly settle into a regime of stellar/anti stellar circulation 
where high altitude winds blow from the dayside to nightside in 
a symmetric pattern with respect to the substellar point (upper 
right panel of Fig. 3; similar to the slowly rotating case of Merlis 
& Schneider 2010). 

To understand these distinct circulation regimes, one needs 
to understand the mechanism responsible for the onset of super- 
rotation on tidally locked exoplanets. Showman & Polvani 



(20 11) showed that the strong longitudinal variations in radiative 
heating causes the formation of standing, planetary-scale equa- 
torial Rossby and Kelvin waves. By their interaction with the 
mean flow, these waves transport eastward momentum from the 
mid latitudes toward the equator, creating and maintaining an 
equatorial jet. As a result, they predict that the strength of the 
jet increases with the day/night temperature difference and that 
the latitudinal width of this jet scales as the equatorial Rossby 
deformation radius 




(6) 



where is the Brunt-Vaisala frequency, H is the pressure scale 
height of the atmosphere (A^ H being the typical speed of gravity 
waves), and /3 = df fdy is the latitudinal derivative of the Corio- 
lis parameter, / s 2 O sin 0, with Q the planet rotation rate and 
9 the latitude (fi - 2Q./Rp at the equator). Roughly, the equa- 
torial Rossby deformation radius is the typical length scale over 
which an internal (gravity) wave created at the equator is signif- 
icantly affected by the Coriolis force. This allows us to define a 
dimensionless equatorial Rossby deformation length 



£ = 



NH 



(7) 



This latter prediction is directly relevant to our case. Indeed, 
although they have rather similar sizes and gravity, G1581 c ro- 
tates around its star in ~ 13 Earth days (d) whereas HD 855 12 b, 
which orbits a brighter K star, has a ~ 58 d orbital period. If 
spin orbit synchronization has occurred, the Coriolis force is 
thus expected to have a much weaker effect on the circulation 
of HD 85512 b. To quantify this, let us evaluate X. in our case. 
In a dry stably stratified atmosphere with a lapse rate dT/dz, the 
Brunt-Vaisala frequency is given by 



r\c„ dz 



(8) 



Cp being the heat capacity of the gas at constant pressure. The 
pressure scale height is given by 



H = 



(9) 



where A;b is the Boltzmann constant and is the mean molec- 
ular weight (in kg) of the air. Assuming an isothermal atmo- 
sphere, this yields 



7^1/2 



(10) 



Assuming a N2 dominated atmosphere, using the equilibrium 
temperature as a proxy for T and disregarding the latitude fac- 
tor, one finds 1. = 1.1 for G1581 c and 2.5 for HD 85512b. In 
our slowly rotating case, the planet is thus too small to contain 
a planetary Rossby wave and the mechanism described above is 
too weak to maintain super-rotation. This also explains why the 
jet is planet- wide in the more rapidly rotating case of G1581 c. 
To confirm this analysis, we have run a simulation of Gl 58 1 c 
where we have, everything else being equal, reduced the rota- 
tion rate by a factor of five so that the X is the same than for 
HD 85512b {Q. - Oorb/5 case; middle column of Fig. 3). As 
expected from the dimensional analysis, in this slowly rotating 
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Fig. 3. Wind maps (arrows) at the 200 mb pressure level (top panel) and near the surface (900 mb pressure level; bottom panel) for G1581 c, 
G1581 c with a slower rotation (see text), and HD 85512 b (from left to right respectively). The mean surface pressure is 1 b for all the models. 
Color scale (same for all panels) and numbers show the wind speed in m/s. Air is converging toward the substellar point at low altitudes in all 
models. However, whereas the standard Gl 58 1 c case clearly exhibits eastward jets (super-rotation) at higher altitudes, both slowly rotating cases 
tend to have stellar/anti stellar winds at the same level. 




Fig. 4. Zonally averaged zonal winds (in m/s) for the dry Gl 58 1 c runs for a total surface pressure of 0. 1 , 1 and 10 bars (from left to right). 



case, no jets appear and we recover a stellar/antistellar circula- 
tion. Another argument in favor of this analysis is the fact the 
location and shape of the cold regions and wind vortices present 
in the G1581 c case (see Figs 1 and 3) are strongly reminiscent 
of Rossby-wave gyres seen in the Matsuno-Gill standing wave 
pattern (Matsuno 1966; Gill 1980; Showman & Polvani 2011). 

Another way to look at this is to consider a multi-way force 
balance (Showman et al. 2013). A way to tilt eddy velocities 
northwest-southeast in the northern hemisphere and southwest- 
northeast in the southern hemisphere, driving super-rotation, 
is to have a force balance between pressure-gradient (V;,(5 = 
A:BArAln p/im^Rp), where AT" is the horizontal day-night tem- 
perature difference and A In p is logarithmic pressure range over 



which this temperature difference is maintained), Coriolis (~ f u 
where u is the wind speed), and advection (^ u^/R^) accelera- 
tion (see Showman et al. (2013) for a discussion on the effect of 
drag). In this case one should have 

V^O~/m~— , (11) 
Rp 

which leads to m ~ fRp and eventually / ~ ^JVp<S>/Rp. 

If / <s; ^JVp<^>/Rp this three-way force balance should re- 
duce to a two way force balance between advection and pressure- 
gradient accelerations which leads to symmetric day-night flow. 
A numerical estimate with the same parameters as above with 
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Ar = 100 K and Alnp = 1 yields a critical rotation period of 
about lOd, close to G1581 c rotation period but much smaller 
than the one of HD 85512 b. This confirms that HD 85512 b 
should have a much more pronounced day-night circulation pat- 
tern than Gl 58 1 c, but it also tells us that Coriolis force is already 
quite weak in the force balance of Gl 581 c explaining the signif- 
icant meridional wind speeds seen in Fig. 3. 

There are a few general differences with previous close-in gi- 
ant planet 3D simulations, however. First, in all our simulations, 
low altitude winds converge toward the substellar point. This re- 
sults from the strong depression caused there by a net mass flux 
horizontal divergence high in the atmosphere driving a large- 
scale upward motion. Second, even at levels where eastward 
winds are the strongest (in the simulations with super-rotation), 
the wind speeds exhibit a minimum (possibly negative) along the 
equator, westward from the substellar point. At a fixed altitude, 
this results from a pressure maximum occurring there because of 
the higher temperature, and consequently higher pressure scale 
height. 
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Fig. 5. Redistribution efficiency rj (ratio of night to dayside outgo- 
ing thermal emission) as a function of atmospheric total pressure (ps) 
for G1581 c (black circles) and HD 85512 b (gray squares), zy = im- 
plies no energy redistribution at all whereas an homogeneous thermal 
emission would yield 77 = 1. The star symbol shows the redistribution 
efficiency of the f2 = f2orb/5 G1581 c case. 

Another feature visible in Fig. 4 is the dependence of the 
strength of the zonal jet to the atmospheric mean surface pres- 
sure (lower pressure causing faster winds). A detailed analysis 
is out of the scope of this study. Nevertheless, these simulations 
seem to confirm the fact that zonal wind speeds increase with 
forcing amplitude (Showman & Polvani 2011). Indeed, an in- 
creasing pressure entails a decreasing thermal contrasts between 
the two hemispheres and a smaller forcing, hence the lower wind 
speeds. 

Surprisingly, the presence of a strong equatorial jet lowers 
the heat redistribution efficiency by the atmosphere. By focus- 
ing the flow around the equator, the Coriolis force prevents an 
efficient redistribution to the midlatitudes on the night side. To 
quantify in a simple way this redistribution, we computed the ra- 
tio 7] of the mean outgoing thermal flux that is emitted by the 
night side over the mean outgoing thermal flux from the day 
side (see Fig. 5). With this definition, no redistribution should 



yield 77 = whereas a very homogeneous emission would pro- 
duce 7] - 1 . As expected, the redistribution increases with total 
pressure. At constant pressure, redistribution seems to be con- 
trolled by the equatorial Rossby number and is more important in 
the slowly rotating cases (HD 855 12 b and O - Oorb/5 Gl 58 1 c 
case). 

Further idealized simulations are needed to fully understand 
the mechanisms controlling both the circulation pattern and the 
energy redistribution on slowly rotating, synchronous terrestrial 
planets. 

3.2.3. Thermal profile 

Finally, in Fig. 6, we present typical temperature profiles ob- 
tained from our simulations (here in the 1 bar case, but profiles 
obtained for different pressures are qualitatively similar). For 
comparison. Fig. 6 also shows the profile obtained from the cor- 
responding global ID simulation (dashed curve), and the spa- 
tially averaged profile from our 3D run (solid curve). 

One can directly see that, in the lower troposphere, temper- 
ature decreases with altitude on the day side following an adi- 
abatic profile and increases on the night side. This tempera- 
ture inversion on the night side is caused by the redistribution 
of heat, which is most efficient near 2-4 km above the surface 
where winds are not damped by the turbulent boundary layer. 
The resulting mean temperature profile cannot be reproduced by 
the ID model, which is to hot near the surface and too cold by 
more than 50 K in the upper troposphere. The discrepancy tends 
to decrease in the upper stratosphere. 
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Fig. 6. Temperature pro- 
files of the 1 bar case for 
G1581c (Grey curves). 
The solid black curve 
stands for the spatially av- 
eraged profile of the 3D 
simulation. For compar- 
ison, the dashed curve 
represents the temperature 
profile obtained with the 
global ID model. 



4. Bistable climate: runaway vs cold trapping. 

As inferred from ID models (e.g. Kasting et al. 1993; Selsis 
et al. 2007) and confirmed by our 3D simulations (not shown), if 
water were present in sufficient amount everywhere on the sur- 
face (aquaplanet case), both our prototype planets (Gl 581c and 
HD 85512 b) would enter a runaway greenhouse state (Kasting 
et al. 1984). The large flux that they receive thus prevents the 
existence of a deep global ocean at their surface. 

However, as visible in Fig. 1, because of the rather ineffi- 
cient redistribution, a large fraction of the surface present mod- 
erate temperatures for which either liquid water or ice would be 
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thermodynamically stable at the surface. In the 10b case, the 
whole night side exhibits mean temperatures that are similar to 
the annual mean in Alaska. It is thus very tempting to declare 
these planets habitable because there always exist a fraction of 
the surface where liquid water could flow. While tempting, one 
should not jump to this conclusion because of the two following 
points. 

• Even if water is thermodynamically stable, it is always evap- 
orating. Consequently, the water vapor is transported by 
the atmosphere until it condenses and precipitates (rains or 
snows) onto the ground in cold regions. This mechanism 
tends to dry warmer regions and to accumulate condensed 
water in cold traps. If these cold traps are permanent, or if 
no mechanism is physically removing the water from them 
- e.g. ice flow due to the gravitational flattening of an ice 
cap, equilibration of sea level (see Sect. 6) - this transport is 
irreversible. 

• Because it is a strong greenhouse gas, water vapor has an 
important positive feedback. The presence of a small amount 
of water can dramatically increase surface temperatures and 
even trigger the runaway greenhouse instability. 

Therefore, even if useful, "dry" simulations cannot be used to 
assess the climate of a planet where water is potentially present, 
and the distribution of the latter^. Hence, in this section we 
present a set of simulations taking into account the full water 
cycle but with a limited water inventory. 

We will show that, as already suggested by Abe et al. (2011), 
climate regimes on heavily irradiated land planets result from 
the competition of the two aforementioned processes and lead 
to the existence of two stable equilibrium states. A "collapsed" 
state where almost all the water is trapped in permanent ice caps 
on the night side or near the poles, and a runaway greenhouse 
state where all the water is in the atmosphere. Furthermore, we 
show that the collapsed state can exist at much higher fluxes than 
inferred by Abe et al. (2011), and that there is no precise flux 
limit, as it can depend on both the atmosphere mass and water 
amount. 



4.1. Synchronous case 

Our numerical experiment is as follows. For our two planets 
and for various pressures, we start from the equilibrium state 
of the corresponding "dry" run described in Sect. 3.2 to which 
a given amount of water vapor is added. To do so, a constant 
mass mixing ratio of water vapor, q^, is prescribed everywhere in 
the atmosphere, and surface pressure is properly rescaled to take 
into account the additional mass. Two typical examples of the 
temporal evolution of the surface temperature and water vapor 
content in these runs is presented in Fig. 7. 

While idealized, this numerical experiment has the advan- 
tage of not introducing additional parameter to characterize the 
initial repartition of water It also provides us with an upper 
limit on the impact of a sudden water vapor release as vapor is 
present in the upper atmosphere where its greenhouse effect is 
maximum. From another point of view, this experiment roughly 
mimics a planet-wide release of water by accretion of a large as- 
teroid or planetesimal and can give us insight into the evolution 
of the water inventory after a sudden impact. 



Because of the sudden addition of water vapor, the atmo- 
sphere is disturbed out of radiative balance. Water vapor sig- 
nificantly reduces the planetary albedo by absorbing a signif- 
icant fraction of the stellar incoming radiation and it also re- 
duces the outgoing thermal flux. The temperature thus rapidly 
increases until a quasi-static thermal equilibrium consistent with 
the amount of water vapor is reached in a few tens of days (see 
Fig. 7). At the same time, water vapor condenses, form clouds, 
and precipitate onto the ground near the coldest regions of the 
surface. As temperature increases, the amount of precipita- 
tion can decrease because the re-evaporation of the falling water 
droplets or ice grains becomes more efficient. The evaporation 
of surface water also becomes more efficient. 
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' At low stellar incoming fluxes, the presence of water also entails a 
positive feedback through the ice albedo feedback so that this argument 
still holds (see e.g. Pierrehumbert 2010) 



Fig. 7. Time evolution of the fraction of water vapor in the atmo- 
sphere over the total water amount and of the mean surface tempera- 
ture for G1581 c with a 200mbar background atmosphere. The black 
solid curve represents a model initialized with a water vapor column 
of 150kg/m^ and the gray dashed one with 250kg/m^. After a short 
transitional period, either the water collapses on the ground in a few 
hundred days (solid curve) or the atmosphere reaches a runaway green- 
house state (dashed curve). 



Then, the system can follow two different paths. If the ini- 
tial amount of water vapor and the associated greenhouse ef- 
fect are large enough, evaporation will overcome precipitation 
and the atmosphere will stay in a stable hot equilibrium state 
where all the water is vaporized (with the exception of a few 
clouds; dashed curve in Fig. 7 and top panel in Fig. 8). Because 
this state is maintained by the greenhouse effect of water vapor 
and that any amount of water added at the surface is eventu- 
ally vaporized, this state is the usual runaway greenhouse state 
(Nakajima et al. 1992). On the contrary, for smaller water va- 
por concentrations, the temperature increase is not large enough 
and the water vapor collapses in the cold traps (solid curve in 
Fig. 7 and bottom panel in Fig. 8). The amount of water vapor 
decreases and the surface slowly cools down again until most of 
the water is present in ice caps in the colder regions of the sur- 
face. An example of the ice surface density of a collapsed state 
at the end of the simulation is presented in Fig. 9. Because subli- 
mation/condensation processes are constantly occurring, this ice 
surface density may continue to evolve, but on a much longer 
time scale that is difficult to model (see Wordsworth et al. 2013 
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Fig. 8. Surface temperature map (in °C) reached at equilibrium for 
the two cases depicted in Fig. 7 (parameters of Gl 58 1 c with a 200 mbar 
background atmosphere). The model shown in the top panel was initial- 
ized with a water vapor column of 250kg/m^ and is in runaway green- 
house state. The bottom panel shows a collapsed state initialized with 
only 150kg/m- of water vapor. 



for a discussion). Then, water vapor is rather well mixed within 
all the atmosphere, and its mixing ratio is mainly controlled by 
the saturation mass mixing ratio just above the warmest region 
of the ice caps. Depending on the simulations, this mass mixing 
ratio varies between 10"'*-10"^. 

To investigate the factor determining the selection of one of 
these two regimes, we have run a large grid of simulations vary- 
ing both the total surface pressure and the initial amount of water 
vapor. As the main effect of water vapor is its impact on the ra- 
diative budget, the quantity of interest is the total water vapor 
column mass (also called water vapor path) 



.={) Jp=0 



(12) 



which is linked to the spectral optical depth of water vapor (rjl) 
by the relation At^ - K^q^pAz - K^Am*. 

Results for the case of G1581 c are summarized in Fig. 10. 
At a given pressure, a small initial amount of water vapor results 
in a collapsed case (blue empty circles). However, there always 
exists a critical initial water vapor path, function of the back- 
ground pressure, above which the runaway greenhouse instabil- 
ity is triggered and develops (red disks). Because redistribution 
efficiency increases with atmospheric mass, cold trapping be- 
comes weaker and the critical initial water vapor path decreases 
with background atmosphere surface pressure. 

We also carried out simulations with a large amount of wa- 
ter (larger than the critical water vapor path discussed above). 
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Fig. 9. Ice surface density of the collapsed state of the 200 mb 
case for G1581 c presented in Fig. 8. The average ice surface density 
is 150 kg/m". Most of the ice has collapsed in the cold traps. 



but distributed only at the surface in the cold traps. For a back- 
ground atmosphere less massive than ~l-5 bar, these simulations 
are found to be stable against runaway. This demonstrates that, 
for a given total amount of water larger than the critical wa- 
ter vapor path, a "moist" bistability exists in the system. Ow- 
ing to the inhomogeneous insolation, the ID notion of critical 
flux above which runaway greenhouse becomes inevitable is no 
longer valid, highlighting the dynamical nature of the runaway 
greenhouse which needs a large enough initial amount of water 
vapor to be triggered. 

When the insolation is decreased, a higher amount of water 
vapor is needed to trigger the runaway greenhouse instability, 
and this at every background pressure. As would be expected, 
planets receiving less flux are more prone to fall in a collapse 
state. However, unexpectedly, results for HD 855 12 b do not fol- 
low this trend. At every background pressure, a lower amount 
of water vapor is needed to trigger the runaway greenhouse in- 
stability, as shown in Fig. 10. This is due to the peculiar fact 
that, because of the higher redistribution efficiency caused by 
the lower rotation rate of the planet, cold traps are warmer (see 
Sect. 3.2 and Fig. 2 for details). As evaporation is mainly con- 
trolled by the saturation pressure of water above the cold traps, 
it is more efficient, destabilizing the climate. 

In addition to the work of Abe et al. (2011) which showed 
that runaway could be deferred to higher incoming fluxes on land 
planets, our simulations demonstrate that, even for a given flux, 
such heavily irradiated land planets can exhibit two stable cli- 
mate in equilibrium. Furthermore, we have shown that the state 
reached by the atmosphere depends not only on the impinging 
flux, but also on the mass of the atmosphere, on the initial wa- 
ter distribution, and on the circulation regime (controlled by the 
forcing, the size of the planet and its rotation rate). This bista- 
bility is due to the dynamical competition between cold trapping 
and greenhouse effect. An important consequence of this bista- 
bility is that the climate present on planets beyond the inner edge 
of the habitable zone will not only depend on their background 
atmosphere's mass and composition, but also on the history of 
the water delivery and escape, as discussed below. 
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Fig. 10. Climate regime reached as function of the initial surface pres- 
sure and water vapor column mass {m*) for G1581 c runs. Red filled 
disks represent simulations that reached a runaway greenhouse states, 
and blue empty circles stand for simulation where water collapsed. The 
blue shaded region roughly depicts the region where collapse is ob- 
served. The gray dashed line represent the demarcation between run- 
away and collapse states for HD 855 12 b. At higher surface pressure, 
redistribution is large and cold traps are less efficient. Less water vapor 
is thus needed to trigger the runaway greenhouse instability. 



4.2. The 3:2 resonance 

Because of the strong tidal interaction with their host star and of 
the small angular momentum contained in the planetary spin, 
rocky planets within or closer than the habitable zone of M 
and K stars are thought to be rapidly synchronized. Indeed, 
standard theory of equilibrium tides (Darwin 1880; Hut 1981; 
Leconte et al. 2010) predicts that a planet should synchronize on 
a timescale equal to 



1 



the synchronization (Makarov & Efroimsky; private communi- 
cation, 2012). Higher order resonances are too weak. 

To assess the impact of a non-synchronized rotation, we per- 
formed the analysis described above in the 3:2 resonance case. 
The rotation rate of the planet is thus three half of the orbital 
mean motion and the solar day lasts two orbital periods. As ex- 
pected, depending on the thermal inertia of the ground and on 
the length of the day, the night side is hotter than in the syn- 
chronized case. More importantly, there is no permanent cold 
trap at low latitudes. However, like for Mercury (Lawrence et al. 
2013), because of the low obliquity, a small area near the poles 
remains cold enough to play this role and ice can be deposited 
there. As the magnitude of the large scale transport of water va- 
por is lower towards the pole compared to the night side (zonal 
winds are stronger than meridional ones), the trapping is less 
efficient. As a result, the transition towards the runaway green- 
house state discussed in the previous section is triggered with a 
lower initial water vapor path. The timescale needed to condense 
all the water vapor near the poles is also much longer. 



5. Evolution of the water inventory 

How stable are the two equilibrium state described above in the 
long term? It is possible that the climate regime exhibited by an 
irradiated land planet changes during its life depending on the 
rate of water delivery and escape of the atmosphere and/or wa- 
ter: triggering runaway when water delivery (or water ice subli- 
mation by impacts) is too large and collapsing water when either 
the background atmosphere or the water vapor mass decreases 
below a certain threshold. Let us briefly put some limits on the 
fluxes involved. 



3 ^GMlRlk2,,At,' 



5. 1 . Atmospheric escape 

When water vapor is present in the upper atmosphere, XUV ra- 
diations can photo dissociate it and the light hydrogen atoms can 
(23-) escape*^. An upper limit on this escape flux is provided by the 
energy limited flux (watson 1981) approximated by 



where rg is the dimensionless gyration radius (r^ - 2/5 for an 
homogeneous interior; Leconte et al. 2011), G is the universal 
gravitational constant, k2,p is the tidal Love number of degree 2 
and Atp is a time lag which characterizes the efficiency of the 
tidal dissipation into the planet's interior (the higher Afp, the 
higher the dissipation). Other variables are defined in Table 1. 
While the magnitude of the tidal dissipation in massive terres- 
trial planets remains largely unconstrained (Hansen 2010; Bol- 
mont et al. 2011), one can have a rough idea of the orders of 
magnitude involved by using the time lag derived for the Eaith, 
fe2,p Afp = 0.305 X 629 s (Neron de Surgy & Laskar 1997). This 
yields a synchronization timescale of 25 000 yr for Gl 581 c and 
10 Myr for HD 855 12 b, preventing these planets from maintain- 
ing an initial obliquity and fast rotation rate (Heller et al. 201 1). 

As demonstrated by Makarov et al. (2012) in the case 
G1581 d, which is more eccentric and further away from its par- 
ent star than G1581 c, tidal synchronization is not the only pos- 
sible spin state attainable by the planet. As for Mercury, if the 
planet started from an initially rapidly rotating state, it could 
have been trapped in multiple spin orbit resonances during its 
quick tidal spin down because of its eccentric orbit. As G1581 c 
also has a non-negligible eccentricity of 0.07±0.06 (Forveille 
et al. 2011), it has almost a 50% probability of being captured 
in 3:2 resonance (depending on parameter used) before reaching 



(kgm ^s '), 



(14) 



where Fxvv is the energetic flux received by the planet and e 
is the heating efficiency. The problem then lies in estimating 
the flux of energetic photons which varies with time (see Sel- 
sis et al. 2007 for a discussion). Estimates of Fg\ are shown 
in Fig. 1 1 where the XUV flux parametrization of Sanz-Forcada 
et al. (2010) and a heating efficiency of e = 0.15 have been used**. 
In some cases, this formula can also be used to estimate the es- 
cape of the background atmosphere itself. 

In atmospheres where the amount of water vapor in the up- 
per atmosphere is limited, the escape flux can be limited by the 
diffusion of water vapor in air. Following Abe et al. (201 1), the 



Of course, the energy available to photo dissociate water molecules 
themselves can also limit escape. Very simple calculations seem to sug- 
gest that this constraint is less stringent in our case, but this conclusion 
strongly depends on the assumed stellar spectral energy density in the 
XUV. 

' To give an idea of the impact on the water inventory, in Fig. 1 1 we 



estimated the lifetime of an Earth ocean by t^, 



where M,,, 



1.4 X 10 kg. This instantaneous estimate does not take 



into account the flux variation during the period considered. 
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Fig. 11. Energy limited escape flux as a function of stellar age 
(e = 0. 15; solid curve) for Gl 581 c (black) and HD 855 12 b (gray). The 
dotted curve is the diffusion limit in the runaway state (/sir(H2) = 0.1 
and Tsu = 300 K), and the dashed curve is an upper limit to the escape 
in the collapse state (/sti(H2) = 10"'* and T^ti = 220 K). The time needed 
to evaporate an Earth ocean of water (a; 1.4 x 10"' kg) is indicated on 
the right scale. 



diffusion-limited escape flux is written as 



^dl = /su(H2)/7i 



(ma - nij) g 



(15) 



where /str(H2) is the mixing ratio of hydrogen in all forms in 
the stratosphere, Tsti is the temperature there, and m-^ and nii 
are the particle mass of air and the species considered, and 
is the binary diffusion coefficient of i in air (for H2 in air bi^- 
1.9xl02i(r/300K)""m-is-'). 

In the collapsed state, the longitudinal and vertical cold traps 
are efficient. Temperature and water vapor mixing ratios are thus 
low in the stratosphere (Tsa ~ 22QK and /str(H2) < lO"'*; dashed 
curve in Fig. 11). On the contrary, in the runaway state, the 
stratosphere temperature is much higher (» 300 K), and water 
is rather well mixed in the atmosphere. The mixing ratio can 
be directly calculated from the total water amount. To show an 
example, we plotted the case of /str(H2) ~ 10"' in Fig. 11. 

Of course, to have a comprehensive view of the problem, 
hydrodynamic and non-thermal (due to stellar winds) escape 
should also be accounted for. However, some conclusions can be 
made. First, it seems that escape of the background atmosphere 
could be efficient during the first few hundred million years to 
one billion year. This comes to support the possibility of a thin 
atmosphere. Second, as expected, escape of water vapor is al- 
ways greater in the runaway state because of the higher tempera- 
ture and humidity in the higher atmosphere (Kasting et al. 1984). 
Third, in the collapse state, cold traps limit the humidity of the 
upper atmosphere and escape is most likely diffiision-limited to 
a rather low value. 

Interestingly enough, the sudden increase of water loss when 
runaway is triggered can provide a strong negative feedback on 
the water content of land planets beyond the classical inner edge 
of the habitable zone. Indeed, at least at early times where the 
XUV flux is strong enough, if a significant fraction of water is 
accreted suddenly and triggers the runaway, water escape will 
act to lower the water vapor amount until collapse ensues. This 
should happen when the atmospheric water vapor content is on 



the order of magnitude of, although lower than, the critical wa- 
ter column amount shown in Fig. 10'". Then, cold-trapped wa- 
ter has a much longer lifetime. This could particularly play a 
role during the early phases of accretion, where a runaway state 
is probable considering the important energy and water vapor 
flux due to the accretion of large planetesimals (Raymond et al. 
2007), a possibly denser atmosphere, and a faster rotation. If the 
planet receives a moderate water amount, the water could be lost 
in a few tens to hundred million years, leaving time to the planet 
to cool down and accrete a significant amount of ice through ac- 
cretion of comets and asteroids later on; progressively or during 
a "late heavy bombardment" type of event for example (Gomes 
et al. 2005). The possibility a such a scenario is supported by 
the fact that a debris disk has been recently imaged in the G1581 
system (Lestrade et al. 2012). 

5.2. Cold trapping rate 

A constraint to the accumulation of ice from such a late water 
delivery is that the delivery rate must be low enough to allow 
trapping to take place without triggering the runaway. To first 
order, our simulations presented in Sect. 4 can be used to quan- 
tify this cold trapping rate. 

Let us consider a case where the water vapor collapses. As 
visible in Fig. 7, except for the very early times of the simula- 
tions, the water vapor content of the atmosphere tends exponen- 
tially towards its (low) equilibrium value. We can thus define an 
effective trapping rate m* - m* /(2 t^) (in kg s^'m"^), where m* 
is the initial water vapor column mass and is the time taken by 
the atmosphere to decrease the water vapor amount by a factor 
of two. This is basically a rough estimate of the critical rate of 
water delivery above which runaway greenhouse would ensue. 




200 400 600 800 

Ps(mbar) 



Fig. 12. Maximum trapping rate per unit surface (m*) reached in our 
simulations as a function of the background surface pressure. The left 
scale shows the total trapping rate at the planetary scale in Earth oceans 
(a 1.4 X 10^' kg) per Myr. 

The maximum trapping rate obtained in our simulations as 
a function of the background atmosphere surface pressure is 



In the numerical experiment shown in Fig. 10, the critical water col- 
umn amount is defined as the minimum water vapor amount to add to 
a cold initial state to force its transition into a hot runaway state. This 
limit would probably be lower if we started with greater temperatures 
has cold trapping would be less efficient. We thus expect that the transi- 
tion from runaway to collapse state occurs at smaller amounts of atmo- 
spheric water vapor 
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given in Fig. 12. As expected, trapping is less efficient when 
the atmospheric mass is higher because cold traps are warmer. 
The exponential decrease can be understood considering that re- 
evaporation efficiency linearly decreases with saturation pres- 
sure that exponentially decreases with temperature following 
Clausius-Clapeyron law. Then, the maximum amount of ice that 
can be accumulated really depends on the duration of the phase 
during which water is delivered to the planet, but the quite high 
trapping rate found in our simulations suggest that if water accre- 
tion occurs steadily, the final water inventory may not be limited 
by the trapping efficiency but by the rate of water delivery after 
the planet formation 

6. On the possibility of liquid water 

Our simulations show that, if the atmosphere is not too thick 
and the rate of water delivery high enough (while still avoiding 
runaway), a significant amount of water can accumulate in the 
cold traps. However, the low temperatures present there prevent 
liquid water stability. Are these land planets possibly habitable? 
Can liquid water flow at their surface for an extended period of 
time? 

A possibility is that both volcanism and meteoritic impacts 
could be present on the night side, creating liquid water, as has 
been recently proposed to explain the evidences of flowing liq- 
uid water during the early Martian climate (Segura et al. 2002; 
Halevy et al. 2007; Wordsworth et al. 2013). Liquid water pro- 
duced that way is however only episodic. 

In this section, we briefly expose various mechanisms 
that could produce long-lived liquid water (as summarized in 
Fig. 13). In particular, we argue that, considering the thickness 
of the ice cap that could be present, gravity driven ice flows sim- 
ilar to those taking place on Earth could transport ice towards 
warm regions where it would melt. Liquid water could thus be 
constantly present near the edge of the ice cap. We also discuss 
the possibility of having subsurface liquid water produced by the 
high pressures reached at the bottom of an ice cap combined to 
a significant geothermal flux. 

We also explored the quasi-snowball scenario proposed by 
Selsis et al. (2007) for Gl 58 1 c where a low equilibrium temper- 
ature is provided by a very high ice albedo (> 0.9) and temper- 
ature slightly exceeds 0°C at the substellar point. Notwithstand- 
ing the difficulty of having such high albedo ices around such a 
red star (Joshi & Haberle 2012), our simulations tend to show 
that this scenario can be discarded. Indeed, the very strong pos- 
itive feedbacks of both water vapor and ice albedo render this 
solution very unstable and no equilibrium state was found. 

6.1. Ice flows 

As described in Sect. 4, if water accumulates on the surface, the 
atmosphere continuously acts to transport it towards the coldest 
regions where it is the most stable (more precipitations and less 
evaporation). As this accumulation seems possible mostly for 
low background atmospheric pressures, there always exists lo- 

' ' Note that this estimation is relevant only for a synchronized planet 
without a thick primordial atmosphere releasing a negligible geother- 
mal flux. All of these assumptions may be far from valid during the 
first few thousand to million years of the planet's life, when the planet 
was still releasing the energy due to its accretion and losing a possibly 
thick primary atmosphere of hydrogen and helium. In addition, the syn- 
chronization timescale estimated in Sect. 4.2 is also of the same order 
of magnitude. 




Fig. 13. Schematic view of the the eyeball planet scenario with a dry 
day side and a thick ice cap on the night side. 



cations at the surface where temperature is below freezing tem- 
perature (see Figs 1 and 2) and only ice is stable. 

If atmospheric transport were the only transport mechanism 
of water on Earth, the latter would also accumulate onto the poles 
and liquid water would only be present when the Sun would melt 
the top of the ice cap during summer. However, other processes 
come into play to limit the thickness of the polar ice caps. Most 
importantly, because ice is not completely rigid, gravity tends to 
drain ice from the upper regions (or regions where the ice cap 
is thicker) towards lower regions creating ice flows, as seen in 
Greenland and Antarctica among other places. These ice flows 
physically remove ice from the cold traps to transport it towards 
regions where it can potentially melt, at the edge of glaciers for 
example. 

Assuming that a sufficient amount of water has been able to 
accumulate on the surface without triggering the runaway green- 
house instability, such ice flows should eventually occur Thus, 
liquid water produced by melting should be present when the 
ice is drained to the dayside where temperatures are above 0°C. 
Quantifying the amount of ice necessary and the average amount 
of liquid water is however much more uncertain, and maybe not 
relevant considering our present (lack of) knowledge concerning 
extrasolar planets surface properties and topography. Let it just 
be said that in greenland, an ice sheet of around 2 km is sufficient 
to transport ice over distances of about 500- 1000 km and create 
ice displacements up to 2-3 km per year near glaciers (Rignot 
et al. 2011; Rignot & Mouginot 2012). As gravity (which is 
the main driver) and typical horizontal length scales both scale 
linearly with radius if density is kept constant, this order of mag- 
nitude should remain valid for more massive land planets. 

As discussed in Sect. 5, ice thicknesses of a few tens of kilo- 
meters could possibly be reached. It is thus possible to imagine 
a stationary climate state where i) a thick enough ice cap pro- 
duces ice flows and liquid water at the edge of warmer regions, 
ii) this liquid water evaporates iii) the resulting water vapor is 
transported back near the coldest regions where it precipitates. 
A question that remains is whether or not the additional amount 
of water vapor that will be released in the atmosphere can be 
large enough to trigger the runaway greenhouse instability that 
would cause the melting and evaporation of the whole ice cap. 

To assess the possibility of this scenario, we have run an- 
other set of simulations with a thick ice cap. Because we do not 
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Fig. 14. Ice distribution for different surface pressure for Gl 58 1 c (ps = 100 mb, 400 mb and 1 b from left to right). Dry regions are in dark gray, 
ice caps are shown in white and liquid water is in blue. For background surface pressures above 300 mb the greenhouse effect of water vapor is 
sufficient to melt water even on the night side. 



want to add any free parameter (namely the total amount of ice) 
and we want a worst-case scenario (the most favorable to the 
triggering of the runaway), we modify our base model so that 
the surface can act as an infinite reservoir of water whose prop- 
erty (liquid or solid) depends solely on the ground temperature 
(above or below 0°C). In particular, we do not take into account 
the latent heat needed to melt water that would have a stabilizing 
effect on the climate by lowering the water vapor evaporation 
rate. 

Another problem comes from the fact that freezing regions 
change as the amount of water vapor increases in the atmo- 
sphere. In our numerical experiment, we thus update the surface 
properties at each time step as follows. A dry surface grid cell 
becomes wet (a reservoir of water) if it has a neighbor having ice 
at its surface: the ice cap extends. As we want to model a case 
where only the edge of the ice cap melts, if a wet surface grid 
cell has no neighbor where water ice is present, it is dried. Then 
the model is initialized from the dry state from Sect. 3.2 where 
a seed of ice has been implanted in the cold trap, and ran until 
either a stationary equilibrium state is reached or the surface is 
completely dried by the runaway greenhouse effect. Test simu- 
lations show that the final state reached does not depend on the 
initial distribution of ice. 

For G1581 c, our simulations show that above 5 bars, run- 
away is inevitable. Below 1 bar, a stationary state is always 
found where ice and liquid water are both present on the surface 
at the same time (see Fig. 14). Depending on the topography near 
the ice cap edge, this could correspond to a situation where lakes 
are constantly replenished by the melting water or where rivulets 
moisten a more extended region. Interestingly, as shown in the 
middle and left panels of Fig. 14, ice does not necessarily cover 
all the nightside because of the radiative effect of water vapor 
which tends to warm the surface, and the inhomogeneous redis- 
tribution of energy discussed in Sect. 3.2. In particular, melting 
first occurs eastward from the dayside near the equator where the 
surface is heated by the hot air coming from the day side through 
the jet, as visible in the 400 mbar case of Fig. 14. 

For HD 85512 b, again, we find that runaway occurs at lower 
background surface pressures. The transition occurs between 
400 and 500 mb. Also, because of the more homogeneous tem- 
perature distribution of the nightside, the transition is not as pro- 
gressive as shown in Fig. 14 for G1581 c. If the 400 mb case 
reaches a steady state with a full ice coverage of the nightside, 
the 500mb case is almost dry. 



These simulations demonstrate that, even in our most strin- 
gent case, in the scenario of a thick ice cap, vigorous ice flows 
and a significant amount of liquid water could be present with- 
out triggering the runaway greenhouse effect. Indeed, in these 
simulations, most of the grid cells next to the ice cap are cov- 
ered by liquid water (see Fig. 14), and the climate remains sta- 
ble. Given our resolution, this corresponds to a wet region with 
a width of about 1000 km, which is probably an overestimation 
of the horizontal extension of the region where melting will oc- 
cur. Interestingly, when there is liquid water at the surface in our 
model, part of the wet surface is always receiving some sunlight, 
because a state with liquid water only on the night side is not sta- 
ble. In this case, water vapor condenses and the water freezes. 
One must bear in mind that, without further information, we did 
not include any runoff of liquid water at the surface, assuming 
that the small amount of flowing water would be evaporated be- 
fore crossing these 1000 km. A more complete model of the ice 
sheet, ice flows and runoff would be necessary to completely set- 
tle this matter 

6.2. Subsurface liquid water? 

Similarly to icy moons such as Europa and Enceladus (Sotin 
et al. 1998; Pappalardo et al. 1999) and to subsurface lakes such 
as the Vostok lake on Eaith (Oswald & Robin 1973), if a thick 
ice cap exists, the possible presence of a significant amount of 
subsurface liquid water at the bottom of the ice layer is also 
worth considering. Indeed, because of the decrease of ice melt- 
ing temperature with pressure and the increase of temperature 
with depth caused by any geothermal flux, it is easier, in a way, 
to melt water below a thick ice sheet than at the surface. 

To assess this possibility, let us estimate the ice thickness 
needed to melt water beneath it as function of the geothermal 
flux. Considering that vertical energy transport in the subsurface 
is mainly performed by thermal conduction, and that quasi-static 
equilibrium is reached, the temperature gradient is given by 



where (f>g is the (positive) geothermal flux and A'th is the thermal 
conductivity of water ice which ranges between 2.22 W/m/K at 
0°C and 3.48 W/m/K at -100°C (citation a trouver). Thanks to 
hydrostatic equilibrium, and assuming a constant density and 
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conductivity for ice (here 920kg/m"' and 0.3W/m/K respec- 
tively), the temperature pressure profile is given by 



Temperature (°C) 

-90. 0. 90. 



Tip) = + 



(p - Ps) , 



(17) 



where and are the local surface temperature and pressure 
and Pice is the mass density of ice. Following Wagner et al. 
(1994), we parametrize the melting curve of water ice by 
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where Tq - 273.16 K and po = 61 1.657 Pa are the temperature 
and pressure of the triple point of water This parametrization is 
valid only between 251.165 and 273.16 K as liquid water is not 
stable below this temperature at any pressure. 

Then, for a given geothermal flux and surface temperature, 
subsurface liquid water can exist if the curves parametrized by 
Eqs. (17) and (18) intersect each other By solving this implicit 
equation, we can find the pressure at which melting occurs and 
define the minimal ice cap thickness needed to reach this pres- 
sure, Zice - (p - Ps)/(piceg) (for the atmosphcrc considered, 
can safely be neglected in the above equations). Numerical es- 
timates of this thickness as a function of surface temperature 
and geothermal flux are shown in Fig. 16 for a typical gravity 
of 18m/s2. 

One can see that even if no geothermal flux is present, sub- 
surface liquid water could exist for warm enough regions. The 
thickness needed would however be rather large (~3-10km) 
which seems unlikely. Depending on the amount of internal flux, 
however, the minimal thickness can be much lower. On Earth, 
radioactive decay and release of primordial heat generate a glob- 
ally averaged geothermal flux of 92 mW/m^ (Davies & Davies 
2010). If we assume that the rate of volumetric heating is con- 
stant, this mean flux should scale linearly with the planet's ra- 
dius, and our more massive planet should produce between 1 .6 
and 1.8 times this flux. With such a flux, the thickness needed 
decreases below 1 .5 km in the coldest regions but could be as 
low as 300 m. 

Finally, another important source of internal flux for close 
in planets is the friction created by the bodily tides raised by 
the host star. An estimate of this tidal dissipation in the syn- 
chronous case is given by Eq. (14) of Leconte et al. (2010) (see 
also Hut 1981). This heating is a function of the orbital ec- 
centricity of the planet and of the efficiency of the tidal dissi- 
pation into the planet's interior (the imaginary part of the Love 
number, k2,pAtp; see Sect. 4.2). For an Earth-like dissipation 
and eccentricities shown in Table 1, the flux can be estimated to 
be ~ 600 mW/m^ for Gl 58 1 c and would contribute in a large 
part to the interior energy budget'^. As shown in Fig. 16, the 
thickness needed to produce liquid water would be much lower 
(< 300 m). Because the measured eccentricity is highly uncer- 
tain, we computed the eccentricity needed to produce the in- 
ternal fluxes shown in Fig. 16 (top scale of the figure). How- 
ever, one must bear in mind that the dissipation constant is also 
highly uncertain, and the tidal flux computed here could easily 
change by an order of magnitude. For HD 85512 b the tidal flux 



Note that this estimation differs from the one of Selsis et al. (2007) 
because they used an eccentricity of 0.16 as initially found by Udry 
et al. (2007). 
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Fig. 15. Surface temperature map (in °C) in the 100 mb case in pres- 
ence of a thick ice cap and liquid water near the terminator (see text). 
Dayside temperatures are up to 50 K higher than in the dry case (upper 
left panel of Fig. 1). 
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Fig. 16. Ice cap minimal thickness (in m) needed to reach melting tem- 
perature at the bottom of the ice sheet as a function of the temperature 
at its surface and internal heat flux. The white area represents the re- 
gion where conditions are too cold for liquid water to exist, even at high 
pressure. This internal flux could be caused by release of primordial 
heat content, radioactive decay and tidal heating. The top scale shows 
the eccentricity needed to provide the given flux for Gl 581 c assuming 
Earth tidal response (see text). 



is much lower (0. 1 mW/m^) because the planet is much further 
away from its parent K star. 

Extended reservoirs of subsurface water can thus not be ruled 
out. The ice thickness needed, if realistic internal fluxes are 
present, is much lower than the one that could potentially be 
accumulated as discussed in Sect. 5. 
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Fig. 18. Synthetic spectrum in reflection (left) and thermal emission (right) during secondary transit of a G1581 c like planet with a 200 mb 
atmosphere at lOparsecs. The black solid and the gray dashed lines respectively stand for the dry collapsed and runaway case presented in Fig. 8. 
Thin dashed lines in the emission part correspond to blackbody spectra at the temperatures indicated. In the runaway case water is abundant in the 
atmosphere and H20 absorption bands are noticeable throughout the spectrum (between 5 and 9 fim in the emission spectrum and above 0.7 /im 
in reflection). The shift of the thermal emission toward shorter wavelengths due to the much hotter surface is also striking. 




Fig. 17. Bond albedo 

as a function of water vapor 
column amount for G1581c 
runs with ps = 200 mb (gray 
dashed curve) and p, = 700 mb 



50 100 150 200 250 (black solid curve). The surface 
2 albedo of 0.3 is shown as the 

'"vap(kg/m ) long dashed horizontal line. 

7. Observable signature 

Can observation allow us to distinguish between the different cli- 
mate regimes presented here? To start answering this question, 
let us consider various possible observables. 

If a given target's flux is too small to use spectroscopy, broad 
band photometry (as has been done for Kepler-lQh (Batalha 
et al. 201 1) and CoRoT-1 b (Leger et al. 2009)) can already pro- 
vide constraints. Indeed, as shown in Fig. 17, the presence of 
water vapor in the atmosphere of a planet orbiting a M dwarf 
strongly decreases its Bond albedo because many water spectral 
bands peak at the same wavelengths than the stellar spectrum. 
However, while this method can hint at the presence of a strong 
absorber in the atmosphere of the planet observed - hinting at 
the presence of an atmosphere itself - a strong degeneracy exist 
between the various possible gases and spectroscopy is needed. 
This signature will of course decrease for more massive, bluer 
stars. 

To go further, we have computed synthetic spectra of various 
models depicted above. We have used the tool developed by 
Selsis et al. (2011) that can model the spectral phase curve of a 
planet in an arbitrary geometrical configuration directly from the 



output of the GCM code. Here we focus only on the spectrum 
that could be observed during the secondary transit of a transiting 
analogue of G1581 c at lOparsecs. An example of the spectra 
computed for the 200 mb runaway and collapse cases presented 
in Figs 7 and 8 are shown in Fig. 18. 

Although the intensity of the bolometric signal is weak, on 
can see that spectral differences between a runaway and a col- 
lapse case are striking. In reflection first, on can easily spot the 
absorption bands of water vapor in the runaway case which are 
not present in the dry, collapse case. Indeed, when water vapor 
becomes a significant constituent of the atmosphere, the latter 
becomes thick in the H2O molecular bands even in the near in- 
frared. A large fraction of the incoming stellar radiation is thus 
absorbed in those bands before being reflected by the surface or 
Rayleigh scattering. Secondly, in the emission spectrum, one 
can see that the planet in runaway can only radiate trough the 
H2O window regions around 4 and lOfim. The surface that is 
heated by the strong resulting greenhouse eff'ect thus radiates at 
shorter wavelengths. 

Further observation may thus, if a suflicient accuracy is met 
(see Selsis et al. (2011) for a more detailed discussion), pro- 
vide strong constraints on the climate regime reached on such 
strongly irradiated land planets as the presence of strong H2O 
spectral features should be correlated with shape of the thermal 
emission. 



8. What's next? 

As for Gl 58 1 c or HD 855 1 2 b, the habitability of planets receiv- 
ing more flux than the classical runaway greenhouse has often 
been discarded on the basis of ID models that predict that an 
ocean cannot be stable under these conditions. 

Here, we have demonstrated that because of the inhomo- 
geneity of the impinging flux and of the natural ability of the 
atmosphere to transport condensable species toward colder re- 
gions where it can be captured, presence of surface water cannot 
be ruled out on the sole basis of the mean absorbed stellar flux. 
In fact, because of the strong positive feedbacks of water vapor, 
two stable equilibrium climate regimes exist for the atmosphere, 
a runaway greenhouse state and a state where all the water is 
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trapped on the surface in cold regions. A key finding of this 
study is that the state in which the atmosphere settles depends 
not only on the incoming flux, but on many factors; atmospheric 
mass, water content and initial distribution, spin state, and circu- 
lation regime. 

If it is not easy to set clear limits on the domain of existence 
of this bistable regime, however, a few trends emerge. First, the 
mechanism described above strongly relies on the existence of 
an eflicient permanent (or slowly evolving) cold trap. The trap- 
ping in a spin orbit resonance or a very low obliquity thus seem 
to strongly favor the onset of the bistability. Second, except for 
changes in the circulation pattern, the domain of existence of 
the collapse state narrows down in the pressure/initial water va- 
por amount parameter space when the incoming flux increases. 
These two argument suggest that the domain of the bistability 
should become smaller when the star becomes more massive 
because regions where tidal interactions causing resonances to 
occur and damping obliquity move at higher fluxes. Again, it 
is difficult to set a clear limit in terms of a spectral type above 
which this bistability cannot occur Let us just say that our sim- 
ulations support the existence of a bistability around M and K 
dwarfs, and that, if not for its thick atmosphere, Venus climate 
would be very close to the one of HD 85512 b, extending this 
statement to G dwarfs. The presence of ice on the poles of Mer- 
cury (Lawrence et al. 2013) also bring observational support to 
this scenario, at least in the limit of a vanishing atmosphere. The 
validity of our results thus seems to extend to a broad range of 
situations. 

Interestingly, because all the parameters affecting the occur- 
rence of a bistable state change during the life of the planet, the 
present climate of a given object may depend on the evolution- 
ary path that it followed. In particular, the history of the water 
delivery and of the atmospheric escape may play a crucial role 
in determining if these objects are dry rocks surrounded by a hot 
steam atmosphere or harbor vast ice deposits on the night side 
and near the poles. 

Another finding is that the atmospheric circulation present 
on synchronous extrasolar planets does not necessary settles in 
a state of super-rotation. As the super-rotation seems to be pow- 
ered by planetary scale Rossby and Kelvin waves, we find that 
it weakens when the wavelength of these waves, i.e. the equato- 
rial Rossby deformation radius, significantly exceeds the plane- 
tary radius. For small and/or very slowly rotating planets (syn- 
chronously rotating on a wide orbit for example) the circulation 
seems to be dominated by a symmetrical stellar to anti stellar 
circulation at high altitudes. 

To go further ahead, observational data are needed. For- 
tunately, exoplanet-characterizing observatories such as EChO 
(Tinetti et al. 2012) and JWST should soon come in line. And 
even if the characterization of Earth analogs around M dwarfs 
will prove very challenging, we will learn a lot from the obser- 
vation of planets beyond the runaway greenhouse limit such as 
transiting versions of the ones described here and hotter Com- 
bining spectral observations of the primary and secondary tran- 
sits (and possibly phase curves) of a number of strongly irradi- 
ated terrestrial exoplanets should indeed help us to unravel the 
variety of atmosphere masses, compositions and histories; so 
far the least constrained parameters used in exoplanetary atmo- 
spheres modeling. 
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